function out = NVtest(ProductionDataNonCum)


out = struct() ;
studlist     = unique(ProductionDataNonCum.sid) ; 
studlisttemp = [] ; %%%%This will be the list of student ids for people who completed >=2 tasks
marginallist = [] ; %%%%This will be the list of student ids for people who completed 1 task but not 2
for ii=1:length(studlist) 
    if max(ProductionDataNonCum.k(ProductionDataNonCum.sid==studlist(ii)))==1 
        marginallist = [marginallist; studlist(ii)] ; %#ok 
    else 
        studlisttemp = [studlisttemp; studlist(ii)] ; %#ok 
    end 
end
id       = ProductionDataNonCum.sid ; 
tau      = ProductionDataNonCum.t ; 
q        = ProductionDataNonCum.k ; 


burnin      = 1*(q==1) ;
logtaudiff  = [] ;
logqdiff    = [] ;
notlast     = [] ;
notlastdiff = [] ;
burnindiff = [] ;
for ii=1:length(studlist)
    indx = find(id==studlist(ii)) ; 
    notlast = [notlast; q(indx)~=max(q(indx))] ; %#ok
    if indx>=2  %%%%Note that if length(indx)<2, then logtaudiff, logqdiff, and burnindiff below will not add any rows...
        tautemp     = tau(indx) ;
        qtemp       = q(indx) ;
        notlasttemp = notlast(indx) ;
        logtaudiff  = [logtaudiff; log(tautemp(2:end)./tautemp(1:end-1))] ; %#ok
        logqdiff    = [logqdiff; log(qtemp(2:end)./qtemp(1:end-1))] ; %#ok
        notlastdiff = [notlastdiff; notlasttemp(2:end)-notlasttemp(1:end-1)] ; %#ok
        burnintemp  = burnin(indx) ;
        burnindiff  = [burnindiff; (burnintemp(2:end)-burnintemp(1:end-1))] ; %#ok
    end 
end 


X          = [burnindiff logqdiff notlastdiff] ;
Y          = logtaudiff ;
tempparams = (X'*X)\X'*Y ;
residuals  = Y - X*tempparams ;
Omega      = diag(residuals.^2) ;
V_HC0      = (X'*X)\(X'*Omega*X)/(X'*X) ;
robust_se  = sqrt(diag(V_HC0)) ;
out.leadcoeff    = tempparams(3) ;
out.sigleadcoeff = robust_se(3) ;
out.t            = out.leadcoeff/out.sigleadcoeff ;
out.CI_95        = [(out.leadcoeff - abs(norminv(0.025))*out.sigleadcoeff) (out.leadcoeff + abs(norminv(0.025))*out.sigleadcoeff)] ;
out.Pvalue       = 2*normcdf(-1*abs(out.t)) ;
out.df           = length(X(:,1)) - length(X(1,:)) ;

% out


end